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-'^This  report  describes  computer  programs  and  their  operational  details 
that  calculate  continuous  seasonal  cycles  of  probabilities  for  values 
of  climatic  variables  and  then  synthesize  distributions  that  will 
describe  probabilistic  levels  of  risk. 
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PREFACE 


This  report  is  one  of  a  series  of  "Research  Reports"  published  by  the 
Integrated  Row  Crop  Management  Systems  Research  Unit  at  the  Southern 
Piedmont  Conservation  Research  Center,  Watkinsville,  GA.  The  purpose 
of  these  reports  is  to  provide  a  mechanism  for  technology  transfer  to 
potential  users  of  information  developed  by  the  scientists  in  this 
Research  Unit.  Below  is  a  list  of  the  reports  developed  in  this 
series.  Any  report  may  be  obtained  by  requesting  it  at  the 
Watkinsville  address  (phone  404/769-5631  or  FTS  250-2425). 
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INTRODUCTION 

We  in  agriculture  readily  recognize  how  dependent  we  are  on 
seasonal  patterns  of  rainfall,  temperature,  streamflow,  and  other 
climatological  variables.  Because  of  this  dependency,  we  are  motivated 
to  develop  methodology  useful  in  the  study  and  treatment  of 
climatological  variates  that  impact  agriculture.  We  wish  to  isolate 
and  define  quantitatively  the  seasonal  variation  of  natural  risk  and 
uncertainty  in  order  to  provide  improved  inputs  into  the  planning  and 
implementation  of  seasonally  dependent  activities. 

We  have  developed  a  series  of  papers  under  the  general  theme  of 
“Stochastic  Impacts  on  Farming"  that  treat  the  determination  of  risk 
and  uncertainty  in  relation  to  resource  management  and  agricultural 
production.  In  the  first  paper,  Thomas  and  Snyder  (1986a)  presented  a 
data  transformation  method  that  was  developed  and  tested  against 
selected  climatological  records.  Then,  Thomas  et  al.  (1986)  published 
a  companion  research  report  that  describes  a  computer  program  which 
transforms  stochastic  data.  In  the  second  paper,  Snyder  and  Thomas 
(1986)  presented  a  method  for  computing  seasonally  continuous 
probability  distributions.  The  third  paper  (Thomas  and  Snyder,  1986b) 
used  simulation  techniques  to  develop  seasonal  variation  of  risk.  This 
research  report  gives  the  computer  programs  and  their  operational 
details  used  in  the  studies  of  paper  two  (analysis)  and  paper  three 
(simulation) . 

In  the  analysis  phase,  we  adapt  two-dimensional  sliding  polynomials 
to  form-free  frequency  distributions  having  mathematically  continuous 
cyclic  variability  through  the  annual  march  of  season.  The  season  is 
defined  by  monthly  samples  of  data  consisting  of  one  observation  per 
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month  for  the  period  of  record.  Such  data  may  be  extremes  of  maxima  or 

minima  or  may  be  monthly  averages  or  totals.  Bi-modal  distributions  or 

statistical  outliers  require  no  special  or  subjective  treatment; 

however,  this  method  is  presently  limited  to  samples  containing  only 

positive  skew.  The  monthly  data  are  standardized  by  a  separate 

transformation  for  each  of  the  twelve  calendar  months.  The 

standardized  data  are  smoothed  by  fitting  the  two-dimensional  sliding 

polynomial  surface  by  least  squares.  Smoothing  is  constrained  by 

imposition  of  boundary  conditions  for  extreme  magnitude  and  for 

maintenance  of  full  seasonal  cyclic  continuity.  Thus,  we  develop 

'  \ 

variate  probabilities  that  are  defined  by  the  mathematical  surface 
based  on  magnitude  of  event  and  month  of  the  year. 

During  the  analysis  phase  the  sliding  polynomial  surface 
establishes  the  relationship  between  magnitudes  of  an  event  and  the 
probability  of  occurrence  of  the  event.  In  the  simulation  phase,  the 
numerical  relationship  is  reversed.  The  surface  is  now  used  to 
calculate  event  magnitudes  from  computer-generated  and  linearly 
distributed  probabilities.  The  process  simply  reduces  to  reverse 
interpolation.  In  the  analysis,  the  derivation  of  form-free  seasonally 
continuous  probability  distributions  establishes  the  stochastic 
characteristics  of  the  climatological  variate  that  includes  seasonal 
cycles.  Simulation  converts  these  characteristics  into  utilitarian 
expressions  of  risk  and  uncertainty. 

ANALYSIS 

The  method  of  analysis  consists  of  two  parts.  First, 
transformation  of  the  variate  is  necessary  to  bring  into  alignment  the 
samples  from  different  months  of  the  year,  factoring  out  the  different 


I 

I 


I! 


I 

I 

n 

II 
II 
II 
II 
II 
II 
II 
II 


6 


means  and  standard  deviations  to  produce  nearly  homogeneous  data  that 
can  be  smoothed  by  a  mathematical  surface.  Resulting  from  this  first 
step  are  12  graphs  of  sample  probabilities  versus  an  abstract  variate 
V.  The  reader  is  referred  to  Thomas  and  Snyder  (1986a)  and  Thomas  et 
al.  (1986)  for  a  full  discussion  of  the  variate  transform  and  a 
computer  program  that  performs  the  transformation. 

The  second  part  of  the  methodology  is  to  smooth  simultaneously  the 
12  monthly  sample  probabilities  and  the  month-to-month  transition  in 
order  to  produce  continuous  probability  functions  that  follow  the 
natural  continuity  of  seasonal  transitions.  A  surface  of  smoothing  is 
required  which  may  be  visualized  as  a  set  of  contours  on  a  cylinder  of 
the  seasons.  The  smoothing  space  is  defined  by  the  magnitude  of  events 
in  v-scale  versus  the  time  by  months.  The  height  of  the  surface  across 
this  space  is  probability.  Two-dimensional  sliding  polynomials  produce 
the  continuous  but  piece-wise  form-free  surface.  Boundaries  of  the 
space  in  v-scale  are  given  by  Snyder  and  Thomas  (1986),  while  in  time 
scale  they  have  a  particular  requirement.  Seasonal  continuity  requires 
that  the  smoothing  surface  pass  through  December  at  the  end  of  the  year 
and  join  with  January  at  the  beginning  of  the  year  with  mathematical 
continuity.  We  used  the  cyclically  continuous  form  of  sliding 
polynomials  as  presented  by  Snyder  (1980). 

The  smoothing  surface  in  the  time  vs  v  space  is  found  by  least 
squares  determination  of  the  values  of  a  grid  of  nodes.  Our  nodal 
arrangement  is  shown  as  Fig.  1.  We  place  lines  of  nodes  at  January, 
April,  June,  October,  and  again  at  January.  By  imposing  boundary 
conditions,  the  initial  5x7  grid  is  reduced  to  a  4x5  grid.  First,  all 
nodes  at  v=0  are  set  to  zero,  and  the  nodes  at  v=-l  are  set  equal  to 
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the  nodes  at  v=l  (Thomas  and  Snyder,  1986a).  Five  nodes  in  v-scale 
remain.  The  right  boundary  at  v=4  is  a  free  boundary.  Next,  the  nodes 
of  the  upper  January  are  set  equal  to  the  nodes  of  the  lower  January  to 
produce  the  cyclindrical  surface.  This  leaves  four  nodes  in  the 
seasonal  time  scale.  This  boundary-constrained  4x5  grid  of  nodes  is 
then  evaluated  by  two-dimensional  least-squares  smoothing  of  the  12 
monthly  samples. 

SIMULATION 

The  method  of  simulation  consists  of  two  numerical  steps:  the 
computer-generation  of  linearly  distributed  pseudo-random  numbers  and  a 
transformation  of  these  numbers  to  other,  different  distributions.  The 
pseudo-random  numbers  are  generated  in  a  scale  from  zero  to  100  and 
each  may  be  considered  a  probability.  Each  value  of  probability  has  an 
equal  chance  of  occurrence.  In  the  second  step,  the  numerical 
relationship  as  developed  in  the  stochastic  analysis  is  reversed  as 
described  earlier.  The  computed  mathematical  surface  is  now  used  to 
calculate  event  magnitudes  from  the  computer-generated  and  linearly 
distributed  probabilities.  Through  the  simulation  process,  we 
synthesize  distributions  that  will  describe  probablistic  levels  of 
risk,  and  that  will  detail  the  variation  of  this  risk  in  a  continuous 
fashion  by  months  throughout  the  year. 

The  complete  process  of  seasonal  simulation  of  risk  requires 
consideration  of  four  elements.  The  first  two,  probability-event 
magnitude  relation  and  the  seasonal  variation  of  this  relationship,  are 
quantified  by  the  least-square  determination  of  the  sliding  polynomial 
smoothing  surface.  These  first  two  elements  are  thus  internal  to  the 
stochastic  analysis.  The  other  two  elements  are  external  and  require 
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subjective  operational  decision.  One  of  these  elements  is  the  length 
of  the  synthetic  samples  that  are  to  be  generated.  This  number  could 
be  described  as  the  length  of  the  planning  period.  Thus,  simulation 
establishes  the  levels  of  risk  of  occurrence  of  events  during  the 
planning  period.  The  last  element  is  the  number  of  synthetic  samples 
to  be  generated.  This  number  sets  a  level  of  precision,  or  a  level  of 
confidence,  that  can  be  placed  on  the  statements  of  risk. 

DISCUSSION 

A  rainfall  sample  was  selected  to  demonstrate  the  results  of  the 
analysis  and  simulation.  This  sample  was  a  99-year  record  (1885-1983) 
of  monthly  total  rainfall  measured  at  Athens,  GA'.  The  monthly  values 
of  mean,  standard  deviation,  skew,  range,  and  transformation  parameters 
for  this  sample  are  given  by  Snyder  and  Thomas  (1986). 

The  smoothing  of  this  monthly  historical  data  set  produces  a 
surface  of  probability.  As  shown  in  Fig.  2,  this  surface  can  be 
displayed  by  contours  of  equal  probability  drawn  in  the  season-event 
magnitude  space  of  the  problem.  The  figure  displays  the  strong  natural 
seasonality  of  rainfall  that  is  typical  of  the  region.  These  contours 
are  calculated  on  the  full  set  of  24-nodes  defining  the  two-dimensional 
surface  as  in  Fig.  1  of  Snyder  et  al .  (1984).  This  includes  the  four 
boundary  nodes  set  to  zero  at  v=0.  The  values  of  these  nodes  in 
probability  scale  and  their  standard  deviations  are  given  in  Table  3  of 
Snyder  and  Thomas  (1986).  The  standard  deviations  are  expressions  of 
localized  variance  (Snyder  et  al . ,  1984)  and  are  alternatives  to  the 
cumbersome  errors  of  prediction  of  conventional  regression  analysis. 
They  provide  measures  of  the  instability  of  the  surface  at  these  nodes 
and  are  one  measure  of  the  goodness  of  fit.  Since  all  deviations  are 
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less  than  1.5%,  we  conclude  that  the  contours  of  probability  are 
closely  defined. 

To  demonstrate  the  simulation  process,  we  chose  the  length  of  the 
planning  period  to  be  20  years.  We  generated  100  samples  for  this 
planning  period.  The  results  of  simulation  are  presented  in  Fig.  3  in 
two  parts.  Part  "a"  shows  simulated  values  of  rainfall  which  is  number 
5  in  ranking  in  the  100  samples  simulated.  Part  "b"  shows  values  which 
is  number  25  in  ranking  in  the  100  samples.  These  figure  parts  are 
thus  simulated  approximations  to  confidence  levels  of  95%  and  75%  since 
risk  values  would  be  exceeded  on  the  average  in  only  5%  or  25%, 
respectively,  of  occurrences  at  selected  risk  levels. 

Fig.  3  shows  event  magnitude  versus  month  for  coordinate  scale. 

The  items  plotted  to  this  scale  are  selected  rank  values  within  the 
planning  period  and  therefore  express  risk.  For  example,  we  see  that 
in  January,  there  is  a  l-in-20  risk  that  monthly  rainfall  will  equal  or 
exceed  about  33  cm.  This  statement  is  made  with  95%  confidence.  These 
contours  of  equal  risk  show  patterns  closely  resembling  the  seasonal 
probabilities  of  Fig.  2.  However,  this  does  not  mean  that  we  should 
perform  the  stochastic  analysis  and  omit  the  simulation.  The 
simulation  is  necessary  to  quantify  risk  and  uncertainty. 

In  this  simulation,  two  elements  of  variation  are  incorporated. 
First,  the  values  of  the  nodes  defining  the  mathematical  surface  are 
varied  in  conformance  with  their  standard  deviations.  After  varying 
the  nodes,  random  numbers  are  generated  and  converted  to  event 
magnitude  using  the  mathematical  surface  through  the  pulsed  nodes. 
Therefore,  both  natural  stochasticity  and  derived  system  variance 
determine  the  variances  of  the  simulations. 
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Fig.  1.  Arrangement  of  nodes  for  smoothing  probabilities  in 
two  dimensions. 
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Fig.  2.  Seasonal  contours  of  probability  for  Athens  monthly 
total  rainfall. 
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Fig.  3.  Seasonal  contours  of  risk  for  Athens  monthly  total 
rainfal 1 . 
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APPENDIX 


The  Appendix  includes  a  description  of  input  and  output  variables, 
program  listing,  and  sample  output  in  the  determination  of  risk  and 
uncertainty  (including  analysis  and  simulation)  of  a  99-year 
monthly-total  rainfall  record. 


The  programs  are  presented  for  the  convenience  of  potential  users. 
While  the  programs  have  been  run  and  tested  on  various  data  sets,  the 
originators  of  the  programs  assume  no  responsibility  for  their  accuracy 
or  adequacy.  Such  responsibility  must  rest  solely  on  the  user.  We 
stand  ready  to  assist  and  advise  within  the  limitations  imposed  by  our 
operating  resources.  The  analysis  program  is  listed  in  Hewlett-Packard 
BASIC  with  Ext.  2.1.  The  simulation  program  is  listed  in  FORTRAN  V. 
(Trade  name  is  included  for  the  benefit  of  the  reader  and  does  not 
imply  an  endorsement  or  preferential  treatment  of  the  named  product.) 
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ANALYSIS  PROGRAM 


This  seasonally  continuous  probability  program  fits  a  two-dimensional 
sliding  polynomial  surface  to  a  sample  made  up  of  one  event  per  month. 
Record  for  each  month  is  converted  to  a  sample  accumulated  probability 
curve  across  40  classes.  Height  of  the  surface  is  probability.  Data  are 
input  from  disk  by  calendar  month  and  classified  in  v-scale  by  month. 


The  program  'WRITEWTINV  provides  data  files  needed  by  the  seasonally 
continuous  probability  program.  The  data  files  cdntain  sliding  polynomial 
weights  and  the  inverse  sums  of  square  matrix.  No  inputs  are  required 
since  it  is  programmed  for  the  fixed  field  of  nodes  as  shown  in  Fig.  1. 
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Analysis  Program  Input  Variables 


Variable  name 
Prtr 

Comment 

Printer  output  device  control. 

(1=CRT,  706=Printer) 

DsnameS 

Sample  input  disk  file  name 

Stnam$ 

Output  disk  file  name 

T5 

Problem  title 

Ny 

Number  of  years  in  record 

Ck$ 

Program  control  variable. 

H(I) 

Monthly  sample  variate 

For  1=1  to  12 

Wt(*) 

Nodal  weights  for  leadt  square  smoothing 

For  *=1  to  41 

B(*,**) 

Inverse  sums-of-squares  matrix  for  least  squares 
For  *=1  to  20:**=1  to  20 

I 

I 

I 

I 

I 

I 

I 

II 

I 

II 
II 
II 


Analysis  Program  Output  Variables 


Variable  name 

Comment 

T$ 

Problem  title 

Hb(I) 

Monthly  sample  mean 

For  1=1  to  12 

Hsd(I) 

Monthly  sample  std.  deviation 

For  1=1  to  12 

Hsk(I) 

Monthly  sample  coeff.  of  skew 

For  1=1  to  12 

Hp 

Minimum  value  of  intermediate  variate,  h 
minimum  class  1 imit 

Cp 

Common  point  of  two  exponential  limbs 
(in  v-scale) 

Vr 

Right  asymptotic  boundary  for. limb  v2  (i 
v-scale) 

F 

Shape  parameter 

D 

Shape  parameter 

Hmin( I ) 

Monthly  minimum  class  limit 

For  1=1  to  12 

Cf(J,I) 

Sample  class  probabilities 

For  1=1  to  41:J=1  to  12 

Lim( J , I ) 

Sample  class  limits 

For  1=1  to  41:J=1  to  12 

Yn(M) 

Sliding  polynomial  solution  nodes 

For  M=1  to  20 

Py(I) 

Smoothed  probability  values 

For  1=1  to  492 

E(I) 

Residual  errors 

For  1=1  to  492 

Nvar(J) 

Nodal  variances 

For  J=1  to  20 

Lu3(J) 

Nodal  sum  of  weights 

For  J=1  to  20 

Ssy 

Total  adjusted  sum  of  squares 

I 
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I 
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II 
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Total  residual  error  square 

Coefficient  of  determination 

Smoothed  probabilities 
For  *=1  to  12:**=1  to  41 
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REM  this  program  FITS  A  2-D  SLIDING  POLYNOMIAL 

REM  #########  SURFACE  TO  A  SAMPLE  MADE  UP  OF  ONE  EVENT  #######?####? 

REM  #########  PER  MONTH.  RECORD  FOR  EACH  MONTH  IS  CONVERTED  ########### 

REM  #########  TO  A  SAMPLE  ACCUMULATED  PROBABILITY  CURVE  ############# 

REM  #########  ACROSS  40  CLASSES.  HEIGHT  OF  THE  SURFACE  IS 

REM  PROBABILITY.  DATA  ARE  INPUT  FROM  DISK  BY  f#fr########## 

REM  #########  CALENDAR  MONTH  AND  CLASSIFIED  IN  V-SCALE  BY 

ncv]  month  #7fiT##rrn-7TTrTrfff f 

REM  #########  CHECK  LEFT  DISK  DRIVE  FOR  NECESSARY  FILES  ############# 

CAT  INTERNAL  4,1" 

LINPUT  "  ARE  ’PROBWTSP'  AND  'PRBINVMATP'  ON  LEFT  DRIVE?  Y 

REM  #########  IF  YES,  CONTINUE  WITH  PROGRAM  ############# 

IF  Ck$="Y"  THEN  200 

REM  #########  IF  NO,  RUN  PROGRAM  TO  GENERATE  FILES  ############# 

DIM  H(IOO) ,Hb(12) ,Hsd(12) ,Hsk(12) ,Lim( 12,41) ,Cf( 12,41)  ,T$[50] ,U( 12) 

DIM  B(20,20),Wt(20),Py(492),E(492),Sy(20),Yn(20),Ab(12,41),Hnnn  12 

DIM  Lul(20) ,Lu2(20) ,Lu3(20) ,Nvar(20) ,Hku( 16) ,Dsname$[30] ,Stnam$[30J 
DIM  Ct(40),Ch(40) 

INTEGER  Ny,I,J,K,N,L,M,Knt 
REAL  Nyl 

INPUT  "TO  PRINT  TO  CRT  ENTER<1>,  HARDCOPY  ENTER<706>" ,Prtr 
LINPUT  "INPUT  SAMPLE  DATA  SET  NAME" ,Dsname$ 

LINPUT  "STORAGE  DATA  SET  NAME",Stnam$ 

CREATE  BDAT  Stnam$,50 
INPUT  "PROBLEM  TITLE", T$ 

OUTPUT  Prtr;T$ 

ASSIGN  OFour  TO  DsnameS 

INPUT  "NUMBER  OF  YEARS  IN  RECORD", Ny 

Nyl=Ny 

LINPUT  "  CAN  HMIN  BE  NEGATIVE  VALUE?  Y  OR  N",Ck$ 

Ks=.2 


OUTPUT  Prtr; 
OUTPUT  Prtr; 
OUTPUT  Prtr; 
OUTPUT  Prtr;" 
REM  ######### 
FOR  1=1  TO  12 
MAT  Ch=  (0) 

MAT  Ct=  (0) 

REM  ######### 
PRINT  I; 

Sh=0. 

Shh=0. 

Shhh=0. 

Hmin(I)=99999. 
FOR  J=1  TO  Ny 


MONTH  MEAN  STD  DEV 

HP  CP  VR  F 

HMIN  " 

BEGIN  MONTH  LOOP 


COEFFS 

SKEW 

D": 


FIND  FIRST  THREE  MOMENTS  FOR  MONTH  I 


############# 


############# 
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PRINT  J; 

ENTER  (aFour;H(J) 

Check=H( J) 

IF  Checl<<Hmin ( I )  THEN 
Hmin( I )=Check 
END  IF 
Sh=Sh+H(J) 

Hh=H(J)*H(J) 

Thh=Thh+Hh 

Th=Th+H(J) 

Shh=Shh+Hh 

Shhh=Shhh+Hh*H(J) 

NEXT  J 
PRINT 

Hb(I)=Sh/Nyl 

M2=(Shh-Sh*Sh/Ny)/Nyl 

Hsd(I)=SQR(M2) 

Hb3=Hb(I)*Hb(I)*Hb(I) 

Hs3=Hsd(I)*Hsd(I)*Hsd(I) 

Hsk(I)=((Shhh/Nyl)-(3*Hb(I)*Shh)/Nyl+(2*Hb3))/Hs3 
REM  #########  KS=0.2  SETS  Hmin  20%  OF  1  STD  DEV  BELOW 
REM  #########  MINIMUM  DATA  VALUE. 

Hmin( I )=Hmin( I )-Ks*Hsd( I ) 

REM  #########  MAY  LIMIT  Hmin  TO  ZERO 
IF  Ck$="N"  THEN 
IF  Hmin(I)<0.  THEN 
Hmin( I)=0. 

END  IF 
END  IF 

OUTPUT  Prtr  USING  850; I ,Hb( I ) ,Hsd( I ) ,Hsk( I ) 

IMAGE  #,5D,7D.4D,2X,6D.4D,2X,3D.4D 
IMAGE  5(5D.4D),6D.2D,2X 

REM  #########  SET  TRANSFORM  PARAMETERS 

Hp=(Hmin(I)-Hb(I))/Hsd(I)+Hsk(I)/8 

Cp=2.+.375*Hsk(I) 

Vr=4.05+(Cp-2.) 

F=L0G((4.-Vr)/(Cp-Vr))/Hp 

D=F*(Vr-Cp)/Cp 

OUTPUT  Prtr  USING  860;Hp,Cp,Vr,F,D,Hmin( I ) 

V1=0. 

REM  #########  CALCULATE  CLASS  LIMITS 

FOR  K=1  TO  39 

Vl=Vl+.l 

IF  Vl>Cp  THEN  1020 

REM  #########  LEFT  SIDE  OF  TRANSFORM 

Ch(K)=-Hsd(I)*(L0G(Vl/Cp)/D+Hsk(I)/8)+Hb(I) 

GOTO  1040 

REM  #########  RIGHT  SIDE  OF  TRANSFORM 

Ch(K)=Hsd(I)*(L0G((Vl-Vr)/(Cp-Vr))/F-Hsk(I)/8)+Hb(I) 
NEXT  K 

Ch(40)=Hmin(I) 

Itt=l 

Cf(I,l)=0. 

Lim(I,l)=99999.99 


############# 

############# 

############# 


############# 


############# 

############# 

############# 
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CALCULATE  SAMPLE  CLASS  PROBABILITIES 


1090  REM  ######### 

1100  FOR  J=1  TO  Ny 
1110  FOR  K=Itt  TO  40 
1120  IF  H(J)>Ch(K)  THEN  1140 
1130  GOTO  1160 
1140  Ct(K)=Ct(K)+l 
1150  GOTO  1170 
1160  NEXT  K 
1170  NEXT  J 
1180  REM 

1190  FOR  J=2  TO  40 
1200  Ct(J)=Ct(J)+Ct(J-l) 

1210  NEXT  J 

1220  FOR  K=Itt  TO  40 
1230  Ct(K)=Ct(K)/Ny 
1240  NEXT  K 

1250  REM  #########  REALIGN  LIMITS  AND  SAMPLE  PROBS  INTO 

1260  FOR  N=1  TO  40 

1270  Cf(I,N-Itt+2)=Ct(N) 

1280  Lim(I,N-Itt+2)=Ch(N) 

1290  NEXT  N 

1300  FOR  N=1  TO  41 
1310  Tp=Tp+Cf(I,N) 

1320  Tp2=Tp2+Cf(I,N)*Cf(I,N) 

1330  NEXT  N 

1340  NEXT  I 

1350  REM  #########  end  OF  MONTH  LOOP 

1360  ASSIGN  OFour  TO  * 

1370  OUTPUT  Prtr;"  SAMPLE  CLASS  PROBABILITIES 

1380  FOR  1=1  TO  12 

1390  OUTPUT  Prtr  USING  1400;I 

1400  IMAGE  #,5X,5D 

1410  NEXT  I 

1420  OUTPUT  Prtr 

1430  IMAGE  #,4D 

1440  FOR  1=1  TO  41 

1450  OUTPUT  Prtr  USING  1430;I-1 

1460  FOR  J=1  TO  12 

1470  OUTPUT  Prtr  USING  1480;Cf(J,I) 

1480  IMAGE  #,3D.3D,3X 
1490  NEXT  J 
1500  OUTPUT  Prtr 
1510  NEXT  I 
1520  OUTPUT  Prtr 

1530  OUTPUT  Prtr;"  CLASS  LIMITS" 

1540  FOR  1=1  TO  12 

1550  OUTPUT  Prtr  USING  1400;! 

1560  NEXT  I 

1570  OUTPUT  Prtr 

1580  FOR  1=1  TO  41 

1590  OUTPUT  Prtr  USING  1430;I-1 

1600  FOR  J=1  TO  12 

1610  IF  Lim(J,I)=99999.99  THEN 

1620  OUTPUT  Prtr;"  ++++++  "; 


41 


CLASSES  #######. 


############# 
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GOTO  1670 
END  IF 

OUTPUT  Prtr  USING  1660 ;Lim( J , I ) 

IMAGE  #,6D.D,2X 
NEXT  J 
OUTPUT  Prtr 
NEXT  I 
OUTPUT  Prtr 

REM  #########  READ  IN  WEIGHTS  A  DATA  POINT  AT  A  TIME  ############# 
REM  #########  ACCUMULATE  XY  VECTOR.  ############# 

ASSIGN  OWrts  TO  "PROBWTSP: INTERNAL ,4, 1" 

FOR  1=1  TO  12 
PRINT  I; 

FOR  J=1  TO  41 
ENTER  @Wrts;Wt(*) 

FOR  L=1  TO  20 
Sy(L)=Sy(L)+Wt(L)*Cf(I,J) 

NEXT  L 
NEXT  J 
NEXT  I 

ASSIGN  OWrts  TO  * 

REM  #########  READ  IN  INVERSE  MATRIX  ############# 

ASSIGN  @Ivn  TO  "PRBINVMATP: INTERNAL, 4,1" 

ENTER  (3Ivn;B(*) 

FOR  L=1  TO  20 
FOR  1=1  TO  20 
Yn(L)=Yn(L)+Sy{I)*B(L,I) 

NEXT  I 
NEXT  L 

REM  ################################################################### 
OUTPUT  Prtr;"SOLUTION  NODES" 

FOR  M=1  TO  8 

OUTPUT  Prtr  USING  1960;M,Yn(M) 

IMAGE  #,4D,2X,3D.5D 
NEXT  M 
OUTPUT  Prtr 
FOR  M=9  TO  16 

OUTPUT  Prtr  USING  1960;M,Yn(M) 

NEXT  M 
OUTPUT  Prtr 
FOR  M=17  TO  20 

OUTPUT  Prtr  USING  1960;M,Yn(M) 

NEXT  M 
OUTPUT  Prtr 
ASSIGN  Olvn  TO  * 

REM  #########  LAY  IN  SLIDING  POLYNOMIAL  SURFACE  THROUGH  ############# 
REM  #########  SOLUTION  NODES.  ############# 

ASSIGN  @Wts  TO  "PROBWTSP: INTERNAL, 4,1" 
lp=0 

FOR  1=1  TO  12 
FOR  J=1  TO  41 
Ip=Ip+l 
PRINT  Ip; 
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ENTER  (3Wts;Wt(*) 

FOR  L=1  TO  20 
Py(Ip)=Py(Ip)+Wt(L)*Yn(L) 

NEXT  L 

REM  #########  CALCULATE  RESIDUALS  ############# 

E(Ip)=Cf(I,J)-Py(Ip) 

Tl=E(Ip)*E(Ip) 

Ses=Ses+Tl 
FOR  L=1  TO  20 
T2=Wt(L)*Wt(L) 

Lu2(L)=Lu2(L)+Tl*T2 

Lul(L)=Lul(L)+T2 

Lu3(L)=Lu3(L)+Wt(L 

NEXT  L 

NEXT  J 

NEXT  I 

ASSIGN  @Wts  TO  * 

OUTPUT  Prtr2340  OUTPUT  Prtr;"  SMOOTHED  PROBABILITY  VALUES 

FOR  1=1  TO  12 

OUTPUT  Prtr  USING  1400;! 

NEXT  I 
OUTPUT  Prtr 
FOR  J=1  TO  41 

OUTPUT  Prtr  USING  1430;J-1 
FOR  I=J  TO  (451+J)  STEP  41 
OUTPUT  Prtr  USING  1480;Py(I) 

NEXT  I 
OUTPUT  Prtr 
NEXT  J 

OUTPUT  Prtr;"  RESIDUAL  ERRORS" 

FOR  J=1  TO  12 

OUTPUT  Prtr  USING  1400;J 

NEXT  J 

OUTPUT  Prtr 

FOR  1=1  TO  41 

OUTPUT  Prtr  USING  1430;I-1 
FOR  J=I  TO  (451+1)  STEP  41 
OUTPUT  Prtr  USING  2550;E(J) 

IMAGE  #,SDD.5D,1X 
NEXT  J 
OUTPUT  Prtr 
NEXT  I 

FOR  L=1  TO  20 
IF  Lu3(L)<=0.  THEN  2630 

Nvar(L)=Lu2(L)/Lul(L)*Ny*12/(Ny*12-16)/Lu3(L) 

GOTO  2640 
Nvar(L)=9999.999 

OUTPUT  Prtr;"  NODAL  VARIANCE";"  NODAL  STANDARD  DEV.";"  SUM  OF  WTS." 

OUTPUT^Prtr^USING  "4D,5D.6D,5D.6D,6D.4D";J,Nvar( J) ,SQR(Nvar(J) ) ,Lu3(J) 
NEXT  J 

Ssy=Tp2-Tp*Tp/492 

OUTPUT  Prtr;"TOTAL  ADJUSTED  SUM  OF  SQUARES  ";Ssy 
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OUTPUT  Prtr;‘'TOTAL  RESIDUAL  ERROR  SQUARE  ";Ses 
OUTPUT  Prtr;"  MEAN  RESIDUAL  ERROR" ;SQR(Ses/492) 


Rs  =  (Ssy-Ses)/Ssy  s 

OUTPUT  Prtr;"CORRELATION  ;Rs;  Lr,/ 

REM  fF########  REORDER  SMOOTHED  VALUES  AND  OUTPUT  TO  DISK 


Knt=0 

FOR  1=1  TO  12 
PRINT  I 
FOR  J=1  TO  41 
Knt=Knt+l 


PRINT  Knt; 

Ab(I,J)=Py(Knt) 

NEXT  J 
PRINT 
NEXT  I 

I 

ASSIGN  @0ut  TO  StnamS 
OUTPUT  @0ut;Cf(*) ,Lim(*) ,Ab(*) 
ASSIGN  @0ut  TO  * 

STOP 

END 


ff############ 
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NODAL  VARIANCE.  NODAL  STANDARD  DEV.  SUM  OF  WTS 

1  .000054  .007335  30.0000 
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"WRITEWTIMV"  Program  Listing 


10 

20 

30 

40 

50 

60 
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90 
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no 
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150 

160 

170 

180 

190 

200 

210 
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240 

250 

260 

270 

280 

290 

300 
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330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 


REM  #####?###################################################### 
REM  ##########  "WRITEWTINV"  ############### 

REM  #########  PROGRAM  TO  WRITE  WEIGHTS  AND  INVERSE  ############# 
REM  #########  SUMS  OF  SQUARES  TO  DISK.  FILES  NEEDED  ############ 
REM  #########  FOR  SEASONALLY  CONTINUOUS  ANALYSIS  ############### 
REM  ##############j?ff##################ff######################### 
OPTION  BASE  1 

DIM  Gc(7,8),C(16),Wt(20),Sx(20,20),Sxinv(20,20) 

INTEGER  Id,Js,Is,L,M,Ic,Iw,Ix,Im,N,K 

REM  ####################j?####################################### 

REM  #########  OPEN  DOUBLE  DO  LOOP  TO  SCAN  ENTIRE  ############### 

REM  #########  12  X  41  SURFACE  MATRIX.  FIND  U  AND  ############### 

REM  #########  Z  VALUES  AND  GET  WEIGHTS  FOR  EACH  ############### 

REM  #########  POINT.  VERTICAL  SCALE  (W)  IS  MONTH  ############### 

REM  #########  HORIZONTAL  SCALE  (X)  IS  V  TRANSFORM  ############## 

REM  ############################################################ 

PRINT  "WEIGHTS  FOR  SEASONALLY  CONTINUOUS  PROB 

PRINT 

Kl=3 

K2=l 

K3=3 

K4=7 

CREATE  BOAT  "PROBWTSP: INTERNAL, 4, 1" ,310 
ASSIGN  OOwt  TO  "PROBWTSP: INTERNAL ,4, 1" 


MONTH  LOOP  BEGINS 
LOCATE  MONTH  ON  SURFACE 


V  TRANSFORM  LOOP  BEGINS 
LOCATE  V  ON  SURFACE 

IS="; 


REM  ######### 

REM  ######### 
ld=0 

FOR  Wm=l  TO  12 
W=l+(Wm-l)/3 
Js=INT(W+. 000001) 

Tw=(W+.000001)-FRACT(W+. 000001) 

U=-.5+W-Tw 
PRINT 

REM  ######### 

REM  ######### 

FOR  Xc=l  TO  41 
PRINT  "JS=";Js;" 

Id=Id+l 

X=l+(Xc-l)/10 

Is=INT(X) 

Tx=X-FRACT(X) 

Z=-.5+X-Tx 
PRINT  Is; 

PRINT  "  ID=  ";Id 

REM  #########  GET  16  WEIGHTS  FROM  SUBROUTINE 
REM  #########  Weight_gen.  PUT  IN  PROPER  4X4 
REM  #########  NODAL  SUB  MATRIX. 

GOSUB  Weight_gen 
FOR  L=1  TO  7 
FOR  M=1  TO  7 
Gc(L,M)=0. 

NEXT  M 
NEXT  L 


############### 

############### 


############### 

############### 


############### 

############### 


I 

I 

I 

I 

fl 

I 

I 


540  lc=0 
550  FOR  L=1  TO  4 
560  FOR  M=1  TO  4 
570  Ic=Ic+l 

580  Gc(L+Js-l,M+Is-l)=C(Ic) 

590  NEXT  M 
600  NEXT  L 

610  REM  #########  SET  BOUNDARIES.  LEFT  BOUNDARY  HAS 
620  REM  #########  ZERO  NODE  AND  ZERO  SLOPE 
630  FOR  Iw=l  TO  7 
640  Gc(Iw,Kl)=Gc(Iw,Kl)+Gc(Iw,K2) 

650  NEXT  Iw 

660  REM  #########  TOP  AND  BOTTOM  BOUNDARIES  WRAP 
670  REM  #########  AROUND  FOR  CYLNDER  FUNCTION. 

680  FOR  Ix=K3  TO  K4 

690  Gc(2,Ix)=Gc(2,Ix)+Gc(6,Ix) 

700  Gc(5,Ix)=Gc(5,Ix)+Gc(l,Ix) 

710  Gc(3,Ix)=Gc(3,Ix)+Gc(7,Ix) 

720  NEXT  Ix 

730  REM  #########  LINEARIZE  THE  Gc(  ,  )  MATRIX,  SAVE 
740  REM  #########  AND  TAKE  THE  SUM  OF  SQUARES. 

750  lm=0 

760  FOR  Iw=2  TO  5 
770  FOR  Ix=K3  TO  K4 
780  Im=Im+l 
790  Wt(Im)=Gc(Iw,Ix) 

800  PRINT  Wt(Im); 

810  NEXT  Ix 
820  NEXT  Iw 
830  PRINT 

840  OUTPUT  @Owt;Wt(*) 

850  PRINT 

860  FOR  N=1  TO  Im 

870  FOR  K=N  TO  Im 

880  Sx(N,K)=Sx(N,K)+Wt(N)*Wt(K) 

890  NEXT  K 

900  NEXT  N 

910  NEXT  Xc 

920  NEXT  Wm 

930  ASSIGN  @Owt  TO  * 

940  REM  #########  SQUARE  THE  TRIANGULAR  SUM-OF- 
950  REM  #########  SQUARES  MATRIX  AND  PRINT 
960  FOR  N=1  TO  Im 
970  FOR  K=N  TO  Im 
980  Sx(K,N)=Sx(N,K) 

990  NEXT  K 
1000  NEXT  N 

1010  PRINT  "  SUMS  OF  SQUARES  MATRIX" 

1020  FOR  N=1  TO  20 
1030  FOR  K=1  TO  10 
1040  PRINT  Sx(N,K); 

1050  NEXT  K 

1060  PRINT 

1070  FOR  K=ll  TO  20 


############### 

#######7^^####### 


############### 

############### 


############### 

############### 


############### 
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1080  PRINT  Sx(N,K); 

1090  NEXT  K 

1100  PRINT 

1110  PRINT 

1120  NEXT  N 

1130  PRINT 

1140  FOR  K=1  TO  8 

1150  PRINT  Sx(K,17); 

1160  NEXT  K 
1170  PRINT 
1180  FOR  K=9  TO  16 
1190  PRINT  Sx(K,17); 

1200  NEXT  K 
1210  PRINT 

1220  REM  #########  USE  SUBROUTINE  Matt_inv  TO  INVERT 
1230  REM  #########  MATRIX. 

1240  GOSUB  Matt_inv 

1250  REM  #########  LOAD  AND  EXECUTE  SEASONALLY 
1260  REM  #########  CONTINUOUS  ANALYSIS  PROGRAM 
1270  PRINT  "  DISK  FILES  HAVE  BEEN  GENERATED  " 

1280  LOAD  "SEAC0NPST",200 
1290  STOP 

1300  Weight_gen:  ! 

1310  IF  U=Us  THEN  1350 
1320  Us=U 
1330  Tr=l 
1340  GOTO  1360 
1350  Tr=2 

1360  ON  Tr  GOTO  1370,1450 
1370  Ul=(((8*U-4)*U-2)*U+l)/256 

1380  U2=(((-32*U+4)*U+8)*U-1)/128 
1390  U3=(((-8*U+4)*U+2)*U-l)/64 

1400  U4=(((32*U-4)*U-8)*U+l)/32 

1410  U5=(((-72*U+36)*U+18)*U-9)/256 
1420  U6=( ( ( 160*U-44)*U-40)*U+11)/128 

1430  U7=(((8*U-4)*U-2)*U+l)/64 

1440  U8=( ( (-96*U+12)*U+24)*U-3)/32 

1450  C(1)=((U4*Z+U3)*Z+U2)*Z+U1 

1460  C(2)=((U8*Z+U7)*Z+U6)*Z+U5 

1470  C(3)=((-U8*Z+U7)*Z-U6)*Z+U5 
1480  C(4)=((-U4*Z+U3)*Z-U2)*Z+U1 

1490  ON  Tr  GOTO  1500,1580 
1500  Vl=(((-24*U+4)*U+22)*U-9)/256 

1510  V2=( ( (96*U-4)*U-40)*U+9)/128 

1520  V3=( ( (24*U-4)*U-22)*U+9)/64 

1530  V4=( ( (-96*U+4)*U+40)*U-9)/32 

1540  V5=(((216*U-36)*U-198)*U+81)/256 

1550  V6=( ( (-480*U+44)*U+296)*U-99)/128 

1560  V7=( ( (-24*U+4)*U+22)*U-9)/64 

1570  V8=(((288*U-12)*U-120)*U+27)/32 

1580  C(5)=((V4*Z+V3)*Z+V2)*Z+V1 

1590  C(6)=((V8*Z+V7)*Z+V6)*Z+V5 

1600  C(7)=((-V8*Z+V7)*Z-V6)*Z+V5 

1610  C(8)=((-V4*Z+V3)*Z-V2)*Z+V1 


############### 

############### 

############### 

############### 
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1620  ON  Tr  GOTO  1630,1710 
1630  Gl=(((24*U+4)*U-22)*U-9)/256 
1640  G2=( ( (-96*U-4)*U+40)*U+9)/128 
1650  63=(((-24*U-4)*U+22)*U+9)/64 

1660  G4=(((96*U+4)*U-40)*U-9)/32 
1670  G5=(((-216*U-36)*U+198)*U+81)/256 
1680  G6=( ( (480*U+44)*U-296)*U-99)/128 
1690  G7=(((24*U+4)*U-22)*U-9)/64 
1700  G8=( ( (-288*U-12)*U+120)*U+27)/32 
1710  C(9)=((G4*Z+G3)*Z+G2)*Z+G1 
1720  C(10)=((G8*Z+G7)*Z+G6)*Z+G5 

1730  C(11)=((-G8*Z+G7)*Z-G6)*Z+G5 
1740  C( 12)={ (-G4*Z+G3)*Z-G2)*Z+G1 
1750  ON  Tr  GOTO  1760,1840 
1760  Hl=(((-8*U-4)*U+2)*U+l)/256 

1770  H2=(((32*U+4)*U-8)*U-1)/128 
1780  H3=(((8*U+4)*U-2)*U-l)/64 
1790  H4=(((-32*U-4)*U+8)*U+l)/32 
1800  H5=(((72*U+36)*U-18)*U-9)/256 
1810  H6=( ( (-160*U-44)*U+40)*U+11)/128 
1820  H7=(((-8*U-4)*U+2)*U+l)/64 
1830  H8=( ( (96*U+12)*U-24)*U-3)/32 
1840  C(13)=((H4*Z+H3)*Z+H2)*Z+H1 

1850  C(14)=((H8*Z+H7)*Z+H6)*Z+H5 

1860  C( 15)=( (-H8*Z+H7)*Z-H6)*Z+H5 
1870  C( 16)=( (-H4*Z+H3)*Z-H2)*Z+H1 
1880  RETURN 
1890  Matt_inv:  ! 

1900  REM  #########  THIS  SUBROUTINE  USES  THE  ADVANCED  ############### 
1910  REM  #########  PROGRAM  FEATURE  TO  INVERT  A  MATRIX  ############### 
1920  CREATE  BOAT  "PRBINVMATP: INTERNAL, 4,1" ,13 
1930  ASSIGN  @Siv  TO  "PRBINVMATP; INTERNAL, 4,1" 

1940  REM  #########  CLEAR  THE  TARGET  MATRIX  ############### 

1950  MAT  Sxinv=  (0) 

1960  REM  #########  AND  PUT  THE  INVERTED  MATRIX  IN  ############### 
1970  MAT  Sxinv=  INV(Sx) 

1980  Check_det=DET 
1990  PRINT 

2000  PRINT  "INVERSE  MATRIX" 

2010  PRINT 
2020  FOR  1=1  TO  20 

2030  FOR  N=1  TO  10 

2040  PRINT  Sxinv(I,N); 

2050  NEXT  N 
2060  PRINT 
2070  FOR  N=ll  TO  20 
2080  PRINT  Sxinv(I,N); 

2090  NEXT  N 
2100  PRINT 
2110  PRINT 
2120  NEXT  I 

2130  REM  #########  CHECK  THE  DETERMINANT  AND  RETURN  ############### 
2140  PRINT  "  DETERMINANT  =  ";Check_det 
2150  BEEP  500,2 
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2160 

WAIT  4 

2170 

OUTPUT 

2180 

ASSIGN 

2190 

RETURN 

2200 

END 

@Siv;Sxinv(*) 
@Siv  TO  * 


SIMULATION  PROGRAM 


This  program  simulates  the  seasonal  variation  of  climatic  risk  by 
randomly  drawing  and  ordering  events  from  a  2-D  sliding  polynomial 
surface.  This  surface  or  smoothing  space  is  defined  by  the  magnitude 
of  events  in  v-scale  versus  time  by  month.  The  height  of  the  surface 
across  this  space  is  probability. 
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Simulation  Program  Input  Variables 


Variable  name 

Comment 

COEF(I) 

Sliding  polynomial  coefficients 
For  1=1  to  4 

AT(I) 

Problem  title 

For  1=1  to  80 

NR,  Nl,  N2 

Random  number  seeds 

NS 

Number  of  samples 

NI 

Planning  period 

HB(I) 

Monthly  mean 

For  1=1  to  12 

HSD{I) 

Monthly  standard  deviation 

For  1=1  to  12 

HSK(I) 

Monthly  coeff.  of  skew 

For  1=1  to  12 

HMIN(I) 

Monthly  minimum  class  limit 

For  1=1  to  12 

YNd.J) 

Solution  nodes 

For  1=2  to  5:  J=3  to  7 

SN(I,J) 

Nodal  standard  deviation 

For  1=2  to  5:  J=3  to  7 
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Simulation  Program  Output  Variables 


Variable  name 

Comment 

AT(I) 

Problem  title 

For  1=1  to  80 

NR,  Nl,  N2 

Random  number  seeds 

NS 

Number  of  samples 

NI 

Planning  period 

HB(I) 

Monthly  mean 

For  1=1  to  12 

HSD(I) 

Monthly  sample  std.  deviation 

For  1=1  to  12 

HSK{I) 

Monthly  sample  coeff.'of  skew 

For  1=1  to  12 

YNd.J) 

Solution  nodes 

For  1=2  to  5;  J=3  to  7 

SN(I,J) 

Nodal  standard  deviation 

For  1=2  to  5:  J=3  to  7 

HP 

Minimum  value  of  intermediate  variate,  h',  at 
minimum  class  limit 

CP 

Common  point  of  two  exponential  limbs 
(in  v-scale) 

VR 

Right  asymptotic  boundary  for  limb  v2 
(in  v-scale) 

F 

Shape  parameter 

D 

Shape  parameter 

H(M,K,N) 

Simulated  variate 

For  N=1  to  12:  M=1  to  NS:  K=1  to  NI 

Simulation  Program  Listing 

C 

c  ############ 
############ 
############ 
############ 
############ 
############ 
############ 


########## 
########## 
########## 
########## 
########## 
########## 
########## 


8110 


8080 


8081 


8010 

1000 

50 

8020 

8022 


SIMULATION  OF  SEASONAL  CONTINUOUS  VARIATION 
OF  CLIMATIC  RISK.  MONTHLY  EVENTS  ARE  DRAWN 
FROM  A  2-D  SLIDING  POLYNOMIAL  SURFACE  AS 
DEFINED  BY  SEASONALLY  CONTINUOUS  PROB. 

ANALYSIS.  HEIGHT  OF  SURFACE  IS  PROBABILITY. 

MONTHLY  EVENTS  ARE  ORDERED  BY  PLANNING 
PERIODS  AND  SAMPLES. 
#################################################################### 
DIMENSION  H(100,50,12),HSK(12),C(4),RN0D(12,7) 

DIMENSION  AT(80),HB(12),HSD(12),HMIN(12) 

DIMENSION  CP(12),HP(12),VR(12),F1(12),D1(12) 

COMMON  YN(7,7),SN(7,7),C0EF(4),SN0D(12,7),YN0D(12,7) 

READ(5,8110)  (C0EF(I),I=1,4) 

F0RMAT(4F12.0) 

READ(5,8080)  (AT( I ) ,1=1 ,80) 

FORMAT (80A1) 

WRITE(6,8081)  {AT( I ) , 1=1 ,80) 

FORMATC  ',80A1) 

READ(5,1000)  NR,N1,N2,NS,NI 
WRITE(6,8010) 

FORMATC  7,20X,'INPUT  DATA'/, 9X, 'SEEDS' ,9X, 'NS' ,9X, 'NI' ) 
F0RMAT(3I5,2I10) 

WRITE(6,50)  NR,N1,N2,NS,NI 
FORMATC  ',6X/,5I8) 

WRITE(6,8020) 

FORMATC  '/,'  MONTH  MEAN  STD  DEV  SKEW  HMIN') 

FORMAT{'  ' ,I4,2F12.3,2F12.3) 

READ(5,1001)  (HB(I),HSD(I),HSK(I),HMIN(I),I=1,12) 

WRITE(6,8022)  ( I ,HB( I ) ,HSD( I ) ,HSK( I ) ,HMIN( I ) , 1=1 ,12) 

KNT=1 


DO  41  1=2,5 
DO  40  J=3,7 

READ(5,8082)  YN( I , J)  ,SN( I ,J) 

8082  F0RMAT(2F10.0) 

YNOD(KNT,J)=YN(I,J) 

SNOD(KNT,J)=SN(I,J) 

40  CONTINUE 
KNT=KNT+3 

WRITE(6,8035)  ( YN( I , J ) ,SN( I , J ) , J=3 , 7) 

8035  FORMATC  ',2F12.3) 

41  CONTINUE 

1001  F0RMAT(4F10.0) 

C  ###########  EXPAND  20  SOLUTION  NODES  TO  12  X  7  #################### 
C  ###########  FIELD  OF  NODES  #################### 

CALL  INTERP 
WRITE(6,332) 

332  FORMATC  '  ,8X , '  HP '  ,8X , '  CP '  ,8X , '  VR '  ,8X , '  F' ,9X , '  D' / ) 

C  ###########  SET  TRANSFORM  ##################################### 

DO  1060  1=1,12 

HP(I)=(HMIN(I)-HB(I))/HSD(I)+HSK(I)/8 

CP(I)=2.+0.375*HSK(I) 

VR(I)=4.05+(CP(I)-2.) 


F1(I)=AL0G({4.-VR(I))/(CP(I)-VR(I)))/HP(I) 
D1(I)=F1(I)*(VR(I)-CP(I))/CP(I) 

WRITE(6,333)  HP( I ) ,CP( I ) , VR( I ) ,F1 ( I ) ,D1 ( I ) 

333  FORMATC  ',5F10.4) 

DO  1060  L=l,7 
RNOD(I,L)=YNOD(I,L) 

1060  CONTINUE 

C  ##########  INITIALIZE  THE  RANDOM  NUMBER  MATRIX 
CALL  RAND(NR,N1,N2,R,1) 

C  ###########  BEGIN  LOOPS: 

C  ###########  M  =  NUMBER  OF  SAMPLE 

C  ###########  I  =  PLANNING  PERIOD 

C  ###########  K  =  MONTH 

DO  1003  M=1,NS 
DO  1004  1  =  1, NI 
DO  1080  K=l,12 
VS=1. 

C  ###########  DRAW  RANDOM  NUMBER  AND  LOCATE  ON  MONTHLY  ############## 
C  ###########  ARC  USING  INTERVAL-HALVING  METHOD  ############## 

CALL  RAND(NR,N1,N2,R,2) 

IF  (R.LT.RN0D(K,3))  GO  TO  1005 

IF  (R.LT.RN0D(K,4))  GO  TO  1006 

IF  (R.LT.RN0D(K,5))  GO  TO  1007 

IF  (R.LT.RN0D(K,6))  GO  TO  1715 

H(M,I,K)=HMIN(K) 

GO  TO  1080 

1005  BL=0.0 
BR=VS 
IS=1 

GO  TO  1008 

1006  BL=VS 
BR=2.0*VS 
IS=2 

GO  TO  1008 

1007  BL=2.0*VS 
BR=3.0*VS 
IS=3 

GO  TO  1008 
1715  BL=3.0*VS 
BR=4.0*VS 
IS=4 

1008  A=(-BL)/(BR-BL)-0.5 

C(l)=(9.*(RN0D(K,IS+l)+RN0D(K,IS+2))-RN0D(K,IS)-RN0D(K,IS+3))/16. 

C(2)=(ll.*(RN0D(K,IS+2)-RN0D(K,IS+l))+RN0D(K,IS)-RN0D(K,IS+3))/8. 

C(3)=(RN0D(K,IS)-RN0D(K,IS+l)-RN0D(K,IS+2)+RN0D(K,IS+3))/4. 

C(4)=(3.*(RN0D(K,IS+l)-RN0D(K,IS+2))-RN0D(K,IS)+RN0D(K,IS+3))/2. 

B=1./(BR-BL) 

1100  Z=A+B*(BL+BR)/2.0 

PC=((C(4)*Z+C(3))*Z+C(2))*Z+C(1) 

IF(R.LT.PC)  GO  TO  1009 

BL=(Z-A)/B 

GO  TO  1010 

1009  BR=(Z-A)/B 

1010  IF(ABS(BR-BL).LT. 0.001)  GO  TO  1011 


################### 

################### 

################### 

################### 

################### 


GO  TO  1100 

1011  V=(BL+BR)/2. 

C  #############  COMPUTE  A  VALUE  H(M,I,K)  #################### 

IF(V.LE.CP(K))  GO  TO  1012 

HS=HSD(K)*(AL0G((V-VR(K))/(CP(K)-VR(K)))/Fl(K)-HSK{K)/8)+HB(K) 

GO  TO  nil 

1012  CONTINUE 

HS=-HSD(K)*(AL0G(V/CP(K))/Dl(K)+HSK(K)/8)+HB(K) 
nil  H(M,I,K)=HS 

C  ###########  PULSE  NODES  WITH  RANDOM  NORMAL  DEVIATE  ################ 
DO  1112  L=3,7 

CALL  RANDN(NR,N1,N2,RANUM) 

RNOD(K,L)=YNOD(K,L)+RANUM*SNOD(K,L) 

1112  CONTINUE 

RN0D(K,1)=RN0D(K,3) 

1080  CONTINUE 

C  ################  END  OF  MONTH  LOOP  #################### 

1004  CONTINUE 

C  ################  END  OF  NI  LOOP  #################### 

DO  1035  K=l,12 
DO  10  J=1,NI 
X6=H(M,J,K) 

DO  20  J1=1,NI 

IF(H(M,J1,K).GE.X6)  GO  TO  20 
H(M,J,K)=H(M,J1,K) 

H(M,J1,K)=X6 

X6=H(M,J,K) 

20  CONTINUE 
10  CONTINUE 
1035  CONTINUE 

C  REMOVE  COMMENTS  ON  CARDS  TO  PRINT  EACH  SAMPLE 

C  WRITE(6,1017)  M 

1017  FORMATC  7, lOX, 'SAMPLE  NO. ',14) 

C  DO  1045  Kl=l,12 

C  WRITE(6,8083)  K1 
8083  FORMATC  '/,10X, 'MONTH' ,14) 

C  WRITE(6,1018)  (K,H(M,K,K1),K=1,NI) 

1018  FORMATC  '  ,8(  I4,F8.3)/) 

1045  CONTINUE 

1003  CONTINUE 

C  ###########  END  OF  NS  LOOP  ##################### 

C  ###########  SORT  SAMPLES  IN  DESCENDING  ORDER  ##################### 
C  ###########  ACROSS  PLANNING  PERIOD  AND  PRINT  ##################### 
DO  2000  01=1,12 
DO  2010  K=1,NI 
DO  2005  1=1, NS 
CHECK=H(I,K,J1) 

DO  2004  0=1, NS 

IF(H(0,K, 01). GE. CHECK)  GO  TO  2004 
H(I,K,01)=H(0,K,01) 

H(0,K,01)=CHECK 

CHECK=H(I,K,01) 

2004  CONTINUE 

2005  CONTINUE 
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2010  CONTINUE 

2000  CONTINUE 

DO  2233  N=l,12 
WRITE(6,2235)  N 

2235  FORMATC  V» lOX MONTH ',14) 

WRITE(6,2223)  (L,L=1,NI) 

2223  FORMATC  SORTED  DATA  ' / ,8X ,20( 12 ,4X) ) 

DO  2001  M=1,NS 

WRITE(6,2222)  M, (H(M,K,N) ,K=1 ,NI ) 

2222  FORMATC  '  ,I3,2X,20F6.2) 

2001  CONTINUE 
2233  CONTINUE 

70  STOP 
END 

C  ###########  SUBROUTINE  RAND  ############### 

C  random  draw  from  a  10  X  10  MATRIX  ############### 

SUBROUTINE  RAND(NR,N1,N2,DRAW,IENT) 

DIMENSION  TAB(10,10) 

IF(IENT.NE.l)  GO  TO  1003 
CALL  RAND0(N1,N12,XRN) 

N1=N12 

II=INT(10.0*XRN)+1 
CALL  RAND0(N2,N22,XRN) 

N2=N22 

JJ=INT(10.0*XRN)+1 
DO  1000  1=1,10 
DO  1001  J=l,10 
CALL  RAND0(NR,NR2,XRN) 

NR=NR2 

1001  TAB(I,J)=XRN 
1000  CONTINUE 


IC=1 

1003  DRAW=TAB(II,JJ) 

CALL  RAND0(NR,NR2,XRN) 

NR=NR2 

TAB(II,JJ)=XRN 

IFC'10D(IC,2).EQ.O)  GO  TO  1005 
CALL  RAND0{N1,N12,XRN) 

N1=N12 

II=INT(10.0*XRN)+1 
GO  TO  1006 

1005  CALL  RAND0(N2,N22,XRN) 

N2=N22 

JJ=INT(10.0*XRN)+1 

1006  IC=IC+1 
RETURN 
END 

SUBROUTINE  RANDO( IX , lY ,YFL) 

IY=IX*65539 

IF(IY)  5,6,6 

5  IY=IY+2147483647+1 

6  YFL=FL0AT(IY)*0.4656613E-9 
RETURN 
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END 

C  #################################################################### 
C  ###########  SUBROUTINE  INTERP  #f?############### 

C  #########frff  EXPAND  20  SOLUTION  NODES  TO  FILL  THE  ################# 

C  ###########  12  X  7  FIELD  OF  NODES  #####?########### 

SUBROUTINE  INTERP 

COMMON  YNN(7,7),SNN(7,7),C0EFN(4),SN0DN(12,7),YN0DN(12,7) 

C  ###ff####  FILL  IN  BOUNDARY  NODES  FOR  CYLINDER  FUNCTION  ############# 
DO  11  J=3,7 
YNN(1,J)=YNN(5,J) 

YNN(6,J)=YNN(2,J) 

SNN(1,J)=SNN(5,J) 

SNN(6,J)=SNN(2,J) 

YNN(7,J)=YNN(3,J) 

SNN(7,J)=SNN(3,J) 

11  CONTINUE 

C  ###########  INTERPOLATE  TO  FILL  IN  GRID  #################### 

KNT=0 

DO  17  1=1,10,3 

DO  16  J=3,7 

YN0DN(I+1,J)=0.0 

YN0DN(I+2,J)=0.0 

SN0DN(I+1,J)=0.0 

SN0DN(I+2,J)=0.0 

VAR1=0.0 

VAR11=0.0 

VAR2=0.0 

DO  15  K=l,4 

YN0DN(I+1,J)=YN0DN(I+1,J)+C0EFN(K)*YNN(K+KNT,J) 

YN0DN(I+2,J)=YN0DN(I+2,J)+C0EFN(K)*YNN{5-K+KNT,J) 

VAR1=VAR1+C0EFN(K)*C0EFN(K)*SNN(K+KNT,J)*SNN(K+KNT,J) 

VAR2=VAR2+C0EFN(K)*C0EFN(K) 

C0EF2=C0EFN(K)*C0EFN(K) 

SN0D2=SNN(5-K+KNT,J)*SNN(5-K+KNT,J) 

VAR11=VAR11+C0EF2*SN0D2 

15  CONTINUE 

SN0DN(I+1,J)=SQRT(VAR1/VAR2) 

SN0DN(I+2,J)=SQRT(VAR11/VAR2) 

16  CONTINUE 
KNT=KNT+1 

17  CONTINUE 

DO  19  1=1,12 
YN0DN(I,1)=YN0DN(I,3) 

SN0DN(I,1)=SN0DN(I,3) 

YN0DN(I,2)=0.0 
SN0DN(I,2)=0.0 
19  CONTINUE 
RETURN 
END 

C  #################################################################### 
C  #########  SUBROUTINE  RANDN  ################ 

C  #########  DRAW  RANDOM  NORMAL  NUMBER  ################ 

SUBROUTINE  RANDN(NR,N1,N2,RANUM) 

RANUM=0. 
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DO  10  1=1,12 
CALL  RAND(NR,N1,N2,X,2) 
RANUM=RANUM+X 
10  CONTINUE 

RANUM=RANUM-6 

RETURN 

END 
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